Quantum crystals in a trapped Rydberg-dressed Bose-Einstein condensate 
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Spontaneously crystalline ground states, called quantum crystals, of a trapped Rydberg-dressed 
Bose-Einstein condensate are numerically investigated. As a result described by a mean- field order 
parameter, such states simultaneously possess crystalline and superfluid properties. A hexagonal 
droplet lattice is observed in a quasi-two-dimensional system when dressing interaction is sufficiently 
strong. Onset of these states is characterized by a drastic drop of the non-classical rotational inertia 
proposed by Leggett [Phys. Rev. Lett. 25, 1543 (1970)]. In addition, an AB stacking bilayer lattice 
can also be attained. Due to an anisotropic interaction possibly induced by an external electric field, 
transition from a hexagonal to a nearly square droplet lattice is also observed. 



Superfluidity which implies a long-range phase coher- 
ence is a crucial property at low temperatures of many 
quantum liquids or gases such as liquid Helium or Bose- 
Einstein condensates (BECs), whereas crystallization im- 
plies a long-range configurational order. Superfluidity 
and crystallization are generally two conflicting proper- 
ties. Penrose and Onsager [1^ were the first to consider 
a BEC in a solid and concluded that such a supersolid 
state simultaneously possessing crystalline and superfluid 
properties was impossible. Since then, this question has 
been revisited by a number of authors [IHl] and has been 
a matter of large speculation for the last forty years. Re- 
cently, the observation of a supersolid phase in ^He sys- 
tems [5] revitalized this fundamental interest. 

An alternative and excellent candidate to study super- 
solidity is in atomic BECs which provide a clean and 
experimentally controllable system. The crystal struc- 
ture in solid helium can be replaced by the modulated 
density in BEC. Density modulated BECs are already 
formed by the imposition of an external potential, cre- 
ating the so-called optical lattices. In these systems, by 
varying the properties of the optical lattice, the conden- 
sate was shown to exhibit a Mott insulator/superfluid 
phase transition [6]-[9j. More recently, it has been shown 
that supersolidity might be present for Rydberg atoms 
in the dipole/van der Waals (vdW) blockade regime [10]- 
fT3] . Cinti et a/. [11 considered a dipole-dipole interaction 
softening at short distance, allowing for a ground-state 
computation that happens to display the properties of 
supersolidity. It proved that a quantum system of inter- 
acting particles can exhibit both crystalline structure and 
superfluidity property. Similar results were obtained by 
Saccani et a/. [12 by using a Heaviside- function interac- 
tion. Based on a mean-field treatment, Henkel et a/. [13 
proposed that a BEC of particles interacting through 
an isotropically repulsive vdW interaction with a soft- 
ened core might support a density modulation. They 
found that the Fourier transform of such interaction has 
a partial attraction in momentum space, which gives 
rise to a transition from a homogeneous BEC to a su- 
persolid phase because of the roton instability (see also 
Refs. [HIIS]). 

Based on the Gross-Pitaevskii (CP) treatment, this 



Letter aims to study the ground states of a trapped 
Rydberg-dressed BEC. Comparing to other CP work 
that did not consider the effect of trapping [13], we 
exactly solve the nonlocal GPE with a trap. In par- 
ticular, we focus on the quasi-2D geometry and the 
strong interacting regime that lead to various fascinat- 
ing ground-state structures of Rydberg-dressed BEC. It 
will be shown that in the supersolid phase, periodic struc- 
tures of Rydberg-dressed BEC can undergo a transition 
from concentric rings to a lattice (or crystalline) if the 
interaction is above some critical value. The lattices are 
formed in terms of crystalline superfluid droplets, called 
quantum crystals, whose onset is characterized by a dras- 
tic drop of the non-classical rotational inertia fraction 
(NCRIF) [4]. The crystalline structure appears to be a 
hexagonal lattice in a quasi-two-dimensional (quasi-2D) 
geometry which can turn into a nearly square lattice if in- 
teraction acquires an anisotropic component in the pres- 
ence of an external electric field (Stark effect). Moreover, 
multilayer crystal structure such as an AB stacking bi- 
layer is also obtained when the frozen axis is relaxed or 
particle number is increased. Crystalline in the case of a 
quasi-one-dimensional (quasi- ID) geometry will also be 
presented. 

We start from a nonlocal Gross-Pitaevskii equation 
(GPE) 
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cally repulsive vdW interaction between Rydberg-dressed 
ground-state atoms [16 with Cq and Rc the effective cou- 
pling constant and blockade radius respectively (we will 
return to the interaction later). Here ^, normalized as 
N = J l^pdr {N is the atom number), denotes the con- 
densate wave function of dressed atoms. M denotes the 
atomic mass and g = Aiih^a/M describes the strength of 
local interaction due to s-wave scattering with scattering 
length a. In the cylindrical coordinates, (p, ^, z), the har- 
monic trapping potential VQ^tir) = [Muj\/2)[p^ + A^z^) 
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FIG. 1. (Color online) Density modulations of the quasi-2D 
Rydberg-dressed condensate with cj = 3, A = 8,7 = 0, and 
a = 10 (a), 50 (b), 2000 (c), and 5000 (d). (e) plots the 
dispersion relation e{kp) [Eq. (|3|] for various interaction cou- 
pling a in a uniform 2D limit. Roton instability occurs when 
ano > 37. (f) shows the calculated NCRIF [Eq. ([5}] with rota- 
tion velocity ujo = 0.01, signifying crystallization at a ^ 2300. 



with radius frequency uj± and aspect ratio A is included 
for simulating real experiments. 

By introducing useful length (Re) and time (r = 
RlM/h) scales, Eq. ^ can be rewritten as 
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where we have redefined the normalized wave function 
t/j = ^/R^/N'^^ the strength of the radius potential 
uj = uj±r^ and the interaction constants 7 = AnNa/Rc 
and a = MNCq/ {h^R^). To obtain ground-state wave 
functions, we computed the governing Eq. ([2| with imagi- 
nary time propagation till the convergence of the normal- 
ized wave function with error less than 10~^. Moreover, 
we have used the method of lines with spatial discretiza- 
tion by the Fourier pseudospectral method. The time 
integration in Eq. ^ is done by the adaptive Runge- 
Kutta method of order 2 and 3 (RK23), which is more 
time efficient due to an adjustable time step. 

For the Rydberg-dressed BEC trapped in a harmonic 
potential, one can define useful characteristic lengths, 
a±_ = ^/h/ {Muj±) and ay = a^/A, corresponding to ra- 
dial and axial potential, respectively. The spectrum and 



the onset of instability are tunable by varying the par- 
ticle number or the confining potential. By varying the 
two ratios Rc/ci± = Rc/ci\\ = >/Aa;, one can effec- 

tively have quasi- ID {^/uo ^ 1 and y/Xuj ^ 1) or quasi-2D 
{^/uj ~ 1 and a/AcJ ^ 1) limits. 

A quasi-2D system with = 3 and A = 8 is first stud- 
ied. Fig. [1] shows the condensate density profile vary- 
ing with the strength of dressing-induced interaction a. 
When a is small, the system displays superfluidity, and 
the ground-state density profile exhibits a central peak 
[see Fig. [TJa)]. As a increases, the central density is too 
high to be stable and thus starts to modulate. Fig. [TJb) 
shows a cratered condensate due to the central insta- 
bility. As a increases further, owing to the roton in- 
stability occurring in the modulated density (see later), 
condensate wrinkles violently and forms a ring structure 
see Fig. [ijc)]. When a is increased above a critical 
value ~ 2300, condensate eventually forms a droplet lat- 
tice. Fig. [TJd) shows the quasi-2D Rydberg-dressed BEC 
forming a hexagonal droplet lattice. 

It is important to check whether the parameter regime 
discussed above is actually experimentally accessible. 
The vdW interaction U{v — v') can be generated from the 
off-resonant dressing of ground-state atoms with high- 
lying Rydberg states. Considering, for example, ground- 
state ^^Rb atoms coupled to excited Rydberg nS state 
^^Rb atoms with n = 60 via a Rabi frequency Vt and a 
red laser detunning A < 0, it will admix a small fraction 
V = (r^/2A)^ of Rydberg character into the ground-state 
atoms for weak dressing <C 1). For two far-distant 
atoms, it leads to an effective interaction Ce/r^ {Cq = 
v'^Cq)^ arising from the strong vdW interaction Cq/v^ 
{Cq = 9.7 X 10^0 a.u.) [n]. At shorter distances, the 
two atoms enter the vdW blockade regime ^ to which 
the effective interaction saturates. Altogether, as given 
earlier, the effective potential between Rydberg-dressed 
ground-state atoms is U{r — r') = CQ/{R^-\-\r — r'\^) with 
the blockade radius Rc = (^6/2^1 A 1)^/^ [13 . As shown, 
the above off-resonant scheme produces only repulsive 
nonlinearities for alkaline atoms [19 . Using typical value 
of Rabi frequency ft = 580KHz and a red detunning 
|A| = 50MHz, blockade radius Rc = 4.5 /im. Moreover, 
the effective lifetime of dressed atoms, l/7eff = l/(^7r)5 
is as large as several seconds with Rydberg state decaying 
rate 7^. ^ 10ms~^ [16 . To achieve the largest coupling 
constant that we are considering, a = 5000, based on 
the above parameters one needs the total atom number 
N = 1.5 X 10^ which corresponds to = vN = 52 
for the number of excited Rydberg atoms. By counting 
the number of droplets ~ 50 in Fig. [ijd), an average of 
one excited Rydberg atom together with 3 x 10^ ground- 
state atoms is within a single droplet. This justifies the 
validity of the CP treatment with the two-body dressing 
interaction /7(r — r') [20] . 

A qualitative understanding of forming the ring-like 
structures [Fig. [ijc)] in a quasi-2D system is given in 
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FIG. 2. (Color online) Comparison of a monolayer (a) and 
a bilayer (b) lattice structures, (a) is the isosurface (isovalue 
0.03) of the density with cj = 3, A = 8, 7 = 0, and a = 5000, 
while (b) is the isosurface (isovalue 0.015) of the density with 
cj = 3, A = 5, 7 = 0, and a = 10000. Inset in (b) shows the 
bilayer being an AB stack. 



the uniform 2D limit. Considering that the system is 
strongly confined in z-direction by a harmonic potential, 
(co'A)^z^/2 but free moving in the xy plane, superfluid 
BEC density can be approximated by n(r) = no(p'^{z) = 
tiq^uoX/t: exp(— cc;A2;^). The 2D excitation spectrum cal- 
culated from the corresponding Bogoliubov-de Gennes 
equations is thus 



e{kp) 
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where kp = |kp| (in units of l/Rc) is mo- 
mentum component in the xy plane and 
U2B{kp) = j U{k)fig{kz)fig{-kz)dkz/2TT with figikz) = 
exp[-i^2/(4^^)] ^ (27rV3)(e-^/Vi^)[e-^/2 - 

2sin(7r/6 - 



/c^) being Fourier 



transforms of the product of (f'g and the scaled interac- 
tion 1/(1 + r^) in Eq. ([2| [22. When deriving Eq. ^ 
and throughout this Letter, contact interaction 7 is set 
to zero. Including 7 will not affect the behaviors of 
roton instability if a is sufficiently large, nor will affect 
the phonon behavior if a is small to which leading term 
of the expansion in U2B{kp) is a positive constant, i.e., 
nonzero 7 only modifies the sound velocity. s{kp) has 
asymptotically a phonon and a free-particle character at 
small and large kp, respectively. However, with U2B{kp) 
having a negative minimum at some finite momentum, 
s{kp) drops near that particular momentum, and eventu- 
ally becomes imaginary when increasing the strength a 
[see Fig. [ije)]. It is estimated that when ariQ > 37, the 
assumed uniform superfluid state is unstable towards 
formation of nonuniform (periodic ring-like) structures. 

To characterize the transition from concentric rings 
to a crystalline hexagonal lattice, we study the non- 
classical rotational inertia fraction (NCRIF), defined by 
(Jo — I)/Io- Here / is moment of inertia of the superfluid 
system under study and Iq is its corresponding classical 
value [5]. As proposed by Leggett [4], NGRIF of the su- 
perfluid system can be calculated under a small rotation. 
In the rotating frame, free energy of a rotating BEG with 
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FIG. 3. (Color online) Comparison of a hexagonal (a) and 
a nearly square (b) lattice structures. Anisotropic ratio k of 
the interaction is 1 (1.4) for (a) [(b)] (see text). 



rotation velocity uq about z axis is 
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where Lz = —i{xdy — ydx) is the z-component angular 
momentum operator and Fq = F{ujq = 0) is the free 
energy of the system without rotation. When uq <^ 1, 
F{ujo) can be expanded as F{ujo) = Fo{iljg) —Iujq/2 with 
t/jg, taken to be real, being the ground state of Fq. Since 
classical moment of inertia is given by /q = / V^^r^dr, we 
obtain for cjo <C 1, 

/ UV -iuJoBz X r) 
NGRIF = 



dr 



uJq J ijj r'^dr 
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where i/j is the ground state of F{ujo). In arriving ([5|, we 
have assumed o:^ i/jg for ujo <^ 1. Therefore NGRIF 
can be obtained by computing Eq. ^ with the solved 

and ipg. Fig. mf) plots the calculated NGRIF as a 
function of the strength a. It is evident that onset of 
crystallization of the BEG droplets is characterized by 
the drastic drop of NGRIF, occurring at a ~ 2300. No 
similar drop appears for the onset of the ring-like su- 
persolid phase occurring at smaller a as a result of full 
rotational symmetry. 

By relaxing the originally frozen z direction potential 
and/or increasing the particle number, density starts to 
modulate in z direction and eventually forms a multilayer 
structure. Figs. |2|a) and [2|b) compare formations of a 
monolayer and a bilayer lattice. The inset in Fig. |2|b) 
indicates clearly that such a bilayer structure is an AB 
stack. 

It is also interesting to note that when Rydberg dress- 
ing interaction becomes anisotropic, hexagonal lattice 
can shift to a nearly square lattice due to distortion 
of the interaction. This can occur when an external 
static electric field E is applied to the system (Stark 
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FIG. 4. (Color online) Density modulations of a quasi- ID 
condensate with cj = 15, A = 0.2, and (a) a = 10, (b) 250, 
(c) 500, and (d) 2000. (e) is the isosurface of modulation (d) 
with isovalue 0.05. 



effect) to which two-photon mechanism will acquire an 
anisotropic component for the interaction, as compared 
to the purely isotropic case with E = Fig. [s] 

shows a nearly (though not perfectly) square lattice 
by using an interaction of the Heaviside-function form 
~ ^(1 — \J ^ K?y'^ + z^) instead of ~ l/(r^ + 1) as 
in Eq. ([2|. Here n corresponds to the anisotropic ratio. 
The static electric field is considered applied along the 
y axis. Without losing the generality, while Heaviside- 
function interaction make the simulation of anisotropy 
more conveniently, it does capture a roton minimum in 
the excitation spectrum. 

The condensate ground state in a quasi- ID system is 
also investigated with = 15 and A = 0.2. Fig. [4] shows 
the modulation of condensate density of a quasi- ID sys- 
tem by varying the strength of dressing-induced interac- 
tion a. When a is small, the system displays superflu- 
idity, and the density has a central peak [see Fig. Qa)]. 
As a increases, the density starts to modulate and has 
multiple peaks. Fig. [4]^b) shows that there are five peaks 
in the condensate, along the axial direction. With a suf- 
ficiently large a, the condensate spontaneously crystal- 
lizes in the axial direction [see Fig. [4]^c)]. As a increases 
further, the condensate starts to modulate in the frozen 
direction, forming a crystalline structure with a central 
hole [see Figs. [4jd) and[4]^e)]. If a is extremely large, 
the condensate starts to cluster in the originally frozen 
direction, and these clusters form a gyroidal chain. 

In summary, based on the Gross-Pitaevskii treatment, 
spontaneously crystalline ground states, called quantum 
crystals, are numerically studied for a trapped Rydberg- 



dressed Bose-Einstein condensate. In a quasi-2D system, 
a hexagonal droplet lattice characterized by a drastic 
drop of the non-classical rotational inertia is shown when 
dressing interaction is sufficiently large. By relaxing the 
originally frozen axis, an AB stacking bilayer lattice is 
observed. We also show that by applying a static electric 
field to make the interaction anisotropic, a nearly square 
droplet lattice can be obtained. 
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